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We present a Monte Carlo study of the finite temperature properties of an extended Hubbard- 
Peierls model describing one dimensional 7r-conjugated polymers. The model incorporates electron- 
phonon and hyperfine interaction and it is solved at the mean field level for half filling. In particular 
we explore the model as a function of the strength of electron-electron and electron-phonon interac- 
tions. At low temperature the system presents a diamagnetic to antiferromagnetic transition as the 
electron-electron interaction strength increases. At the same time by increasing the electron-phonon 
coupling there is a transition from a homogeneous to a Peierls dimerized geometry. As expected 
such a Peierls dimerized phase disappears at finite temperature as a result of thermal vibrations. 
More intriguing is the interplay between the electron-phonon and the electron-electron interactions 
at finite temperature. In particular we demonstrate that for a certain region of the parameter space 
there is a spin-crossover, where the system transits from a low-spin to a high-spin state as the tem- 
perature increases. In close analogy to standard spin-crossover in divalent magnetic molecules such 
a transition is entropy driven. Finally we discuss the role played by the hyperfine interaction over 
the phase diagram. 



I. INTRODUCTION 

Over the past few years the field of organic spin- 
tronics has received interest from a growing and di- 
verse scientific communitjff^. This is primarily driven 
by its potential for opening new avenues to cheap, 
low- weight, mechanically flexible, chemically inert and 
bottom up fabricated spin-devices. The most im- 
portant advantage of organic material^^''^ over their 
inorganic counterparttP^ is the strong suppression of 
any efficient spin flipping mechanisms, resulting in ex- 
tremely long spin lifetimes. Spin orbit coupling is 
very weak in organic semiconductors, as Carbon has 
a low atomic number {Z = 12) and the strength of 
spin orbit interaction is proportional to Z*. Further- 
more also the hyperfine interaction is weak in molecules, 
since there are very few nuclei with non-vanishing mag- 
netic moment in the upper part of the periodic ta- 
ble and the molecular orbitals typically responsible for 
the electron transport are extended 7r-typem. A recent 
experimental workP on (La,Sr)Mn03(LSMO)/ins[8- 
hydroxyquinoline] Aluminum(Alq3)/Co tunnel junctions 
reported a giant tunnel magneto resistance of up to 300% 
at low temperatures. This is a value that compares well 
with the best inorganic tunnel junctions^"^ and it is at- 
tributed to a favorable interaction betwee n the organic 
molecules and the magnetic metallic surfacd^^'^^. 

Amongst the various possible materials for organic 
spintronics, 7r-conjugated organic semiconductors ap- 
pears very appealin g^^-^- '. This is because of their very 
long and relatively temperature independent spin relax- 
ation time and their ability to form good interfaces with 
metal electrodes when incorporated in spin- valve-like de- 
vices. The most relevant structural feature of such poly- 
mers is their planar shape. The 7r-electron wave function 



distends from the molecular plane and it can easily inter- 
act with the wave function of adjacent molecules. There- 
fore a face-to-face molecular configuration is usually sta- 
bilized through strong tt-tt bonding and diffusive Van der 
Waals interaction. This leads to a molecular stacking ar- 
rangement resulting in the formation of a low dimensional 
lattice. Because of this peculiar structural conformation 
many of the 7r-conjugated polymers can be described by 
simplified one dimensional (ID) model HamiltonianJi^. 

An important aspect for the successful integration of 
organic semiconductors in magnetic memories and in 
magnetic switching devices is the feasibility of manip- 
ulating the spin orientation in the organic media. This is 
difficult to achieve in non-magnetic molecules by the tiny 
non-equilibrium spin population originating from spin in- 
jection. In fact the standard techniques for manipulation 
used in inorganic semiconductors, for example optical 
methods, are ineffective because of the weak spin-orbit 
interaction. A more promising option is that of manipu- 
lating the internal spin degrees of freedom of the organic 
medium, when this is magnetic. Intriguingly there is a 
vast class of molecules, generally known as spin crossover 
compounds, whose spin state can be changed from low 
spin to high spin by an external perturbatiorliS. Since 
the crossover transition is entropy driven, it is most typ- 
ically achieved by increasing the temperature, although 
also light, pressure and electro-chemical redox reactions 
can all produce it. Most recently the possibility of spin 
crossover driven by static electric fields has been pro- 
posed theoretically^^ Notably all these processes oc- 
cur in compounds incorporating a transition metal [usu- 
ally Fe(II)], which is responsible for the magnetic mo- 
ment. Thus, it is interesting to explore the possibility 
of obtaining spin crossover in organic 7r-conjugated poly- 
mers. 
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In this work we use a combination of energy minimiza- 
tion techniques and Monte Carlo (MC) simulations to 
investigate the temperature dependent phase diagrams 
for a model Hamiltonian describing 7r-conjugated poly- 
mers. In particular, we explore the parameter space of 
the model in the search of a spin crossover transition. 
Intriguingly we find that the interplay between on-site 
Coulomb repulsion and electron-phonon (el-ph) coupling, 
leading to Peierls distortion, can be responsible for spin 
crossover. 



II. COMPUTATIONAL METHODS 

We consider an extended single-site Hubbard-Peierls 
modeP2for a ID lattice described by the following Hamil- 
tonian 

+ n'i^n^i + Jh X S'i • [c|„((Tq/3)Q;3] + 

where cl^{cia) denotes the creation (annihilation) oper- 
ator for an electron at site i with spin a, tic = e is the 
onsite energy and tij is the transfer integral for an uni- 
form undimerized lattice. Here we consider only nearest 
neighbour hopping, i.e. tij = t for i = j ± 1 and Uj = 
for any other pair. The other microscopic parame- 
ters of the model are the el-ph coupling parameter, a, the 
Hubbard repulsion strength, U, the hyperfine exchange 
Jh, and the elastic constant, k. Thus, the second term 
of the Hamiltonian, in addition to electron hopping, de- 
scribes the el-ph coupling, with qi being the atomic dis- 
placement of site i (we consider ID longitudinal motion 
only). The third term is the standard on-site Hubbard 
repulsion while the fourth one describes a Heisenberg- 
like interaction between the electron spins, cl^{cTai3)cip, 
and a set of classical vectors, 5*^, representing the nuclear 
spins. Finally the last term is classical elastic energy. 

In all our calculations we consider only the half-filling 
case (one electron per site) so that the model exhibits 
an insulating behaviour at T = for an undistorted infi- 
nite chain. The data presented here are for 10 site long 
chains as further tests for longer chains (20 sites) gave 
rather similar results. Our approach consists in replac- 
ing the Hubbard term with its unrestricted mean-field 
approximation^"^ and then in solving the Hamiltonian 
self-consistently for different lattice displacements {qi}- 
Energy minimization is performed by conjugate gradi- 
ent over {qi} and further verified by additional simulated 
annealing^^. The bond lengths, Xi, are calculated from 
the ground state displacements as Xi = d + g^+i — qi, 



where d is the equillibrium bond distance. The main ob- 
servables calculated are the dimerization parameter, D, 
and the local magnetic moment per site, m^. These are 
defined respectively in Eq. ^ and Eq. ([3| 

D ^ ^ \xi - x.i+i \ + \xi+i - Xi+2\ , (2) 

4 = 1 

m,^\nl-nj\. (3) 

Note that in Eq. ([3| nf is the electron spin density at site 
i, which is calculated as the diagonal matrix elements of 

occupied 

the density matrix p = {|*„)(^'„|^3e^}, where 

ri=l 

En and |5'n) are respectively n-th electronic eigenval- 
ues and eigenvectors. Our results are then plotted in a 
phase-diagram like form, where the different phases are 
presented as a function of a and U. 

For finite temperature simulations we consider the sys- 
tem described by the Hamiltonian of Eq. ([T]) and evolve 
the classical dynamical variables {ojl and {Si} by using 
the standard Metropolis algorithrrPSl. Note that since 
Jh ^ i the hyperfine interaction has little effect on the a- 
U phase diagram, so that in what follows we will neglect 
it unless otherwise indicated. As such the only classical 
dynamical variables are the atomic displacements {qi}. 
In the Metropolis algorithm the acceptance probablity of 
a new state is unity if the new configuration has an en- 
ergy lower than that of the old configuration. Otherwise 

it is given by the Boltzmann factor e ''b'^ , where AG 
is the difference in the Gibbs free energy between the 
old and new configuration. By using the grand-canonical 
ensemble the Gibbs free energy G can be written as 

n=l 

where the chemical potential fi is obtained from 

^ " E e^}[£;,.({?.})-A'] + 1 ' 

n 

with (3 = -j^^ being the inverse temperature and N the 
total number of electrons {N = 10 here). For every value 
of {a/t, U/t) and each temperature the system is allowed 
to reach equilibrium. Then both D and mi are calculated 
over several million MC steps. 

In what follows we will express all energy related quan- 
tities (including the temperature) as a function of i, 
which sets the energy scale of the problem. The on-site 
energy is taken to be zero and A; is 5 t/ k? . Note that for 
t ~2.5 eV this corresponds to k = 12.5 eV/A^, which is 
in between the value for the H2 molecule and that of Au 
monoatomic chainJ^. Also it is important to note that 
k > 5t/A^ is a value commonly use d in r ecent literature 
about transport in organic polymer j^^l^^l 
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III. RESULTS AND DISCUSSION 
A. Ground State 

Let us begin our discussion by investigating the a-U 
phase diagram at T = 0, which is presented in Fig. [T] 
This is populated by four difference phases characterized 
by the different combined structural and magnetic prop- 
erties of the chain. In particular there are two magnetic 
states and two geometrical configurations. For small a/t 
and U jt the chains are undistorted (the atomic spacing 
is approximately uniform throughout the chain) and in a 
non-spin-polarized diamagnetic (DM) state. We denote 
this phase as DM-undistorted. As U jt increases eventu- 
ally a spin-polarized solution develops. This is the one 
expected from the mean field Hubbard model at half- 
filling, i.e. it is a antiferromagnetic (AFM) phase, where 
local moments form at each atomic site, but their orien- 
tation alternates along the chain. Such an AFM phase 
may or may not be accompanied by a structural distor- 
tion, depending on the value of ajt. In general however 
there is a vast region of the model parameter space, where 
no significant distortion appears for the AFM spin state. 
We denote this phase as AFM-undistorted. 

As ajt increases for moderate V jt the system progres- 
sively developed a Peierls instability and makes a transi- 
tion to a geometry where long bond distances alternate to 
short ones. This dimerized phase, expected for the non- 
interacting case, remains diamagnetic (DM-dimerized) 
for small U jt but can coexist with a AFM solution for 
a significant range of parameter. In summary the phase 
diagram is characterized by a competition between the 
on-site repulsion, driving the magnetic instability, and 
the el-ph coupling, driving the Peierls distortion. 

In the discussion we have assigned the phase boundary 
of the magnetic transition to the condition 1"^*! 0- 
In contrast assigning the phase boundary to the Peierls 
transition is more complicated since D changes contin- 
uously upon increasing a. Thus we have used the op- 



erational definition of placing the phase boundaries at 
D — 0.17, which is interpreted as representative of strong 
dimerization. The complete evolution of D as a function 
of U jt and a/i is presented in the right hand side panel 
of Fig. [T] The figure clearly reveals the interplay be- 
tween el-ph coupling and Coulomb repulsion. In fact D 
grows almost linearly with ajt for small U jt but then is 
drastically reduced asU/t grows. 



B. Finite temperature phase diagram 

We now move to study the finite temperature prop- 
erties of our model. These are summarized in Fig. [2] 
where we present the a-U phase diagram for two repre- 
sentative temperatures, respectively (5 = 200,400 1/eV 
(these correspond respectively to T = 0.002 0.001 1 or, 
for t = 2.5 eV, to T = 58,29 K). The most important 
feature of the finite temperature plots is the complete ab- 
sence of a structural phase transition. This means that in 
general the system does not dimerize any longer as the 
temperature is increased. The dimerization is instead 
replaced by a general increase in bond length and by a 
random distribution of the various bonds along the chain. 

More details about the structure of the chains at finite 
temperature can be found in the right-hand side panes of 
Fig.[2j which show the quantity D [defined in Eq. ([2|] as a 
function of U/t and a/t for T = O.OOlt and 0.002i. Note 
that at finite temperature the quantity D plotted is a 
measure of the disorder in the bond distances of the chain 
as thermal vibrations onset. As such D is maximum at 
high values of the el-ph coupling and low U/t but falls 
rapidly asU /t increases. 

The second most striking feature of the finite tempera- 
ture a-U phase diagram is the movement of the DM- AFM 
phase boundary towards lower U/t as the temperature 
increases. This essentially means that as the tempera- 
ture grows it takes less on-site Coulomb interaction to 
drive the system towards a magnetic instability. We fur- 
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FIG. 2: (Color online) Finite temperature a-U phase diagram for the model ID polymers investigated. In the left panel we 
present the schematic phase diagram for two different temperatures, respectively T = 0.001 1 and T = 0.002 t. Note that now 
there are only two phases (DM and AFM) and there is no longer a distorted (Peierls) geometrical configuration. In the picture 
we also report the DM-AFM phase boundary for the T = case. Thus the region comprised between the solid black line either 
the solid red of the dashed blue line is characterized by a DM to AFM phase transition (spin crossover - SC - region) as a 
function of temperature. In the middle (T = 0.001 1) and the right-hand side panel (T = 0.002 t) we present the parameter D 
as a function of a/t and U/t which indicates the degree of disorder in the bond distances at finite temperature. 



ther explore this finding in figure [3] wliere we present 
the critical U value, Uq, at which the magnetic phase 
develops. This effectively represents the position of the 
DM-AFM phase boundary. In particular Uq is plotted as 
a function of the temperature, T, and for three different 
el-ph strengths a/t. 




0.001 0.002 0.003 0.004 



FIG. 3: (Color online) Critical on-site Coulomb repulsion, Uc, 
needed for a magnetic solution as a function of temperature 
and for three different values of the el-ph coupling strength. 
The error bars correspond to spacing of the interpolation grid 
used to extract Uc- 



DM-AFM phase boundary as the temperature increases 
is the fact that there is a vast region in the a-U diagram 
where the system makes a DM to AFM transition as the 
temperature increases. Such a region is the one enclosed 
between the two vertical lines marking respectively the 
phase boundary at finite temperature and at T = in 
figure [2j For the particular values of a and U found in 
such a region (called the spin crossover -SC- region) there 
is a temperature driven spin crossover. This is analyzed 
next. 



C. Spin crossover 

We now explain the spin crossover transition by us- 
ing the standard framework of spin crossover usually 
employed for magnetic molecules incorporating divalent 
transition metalglSI. In general the thermodynamically 
stable phase at finite temperature of a system that can 
assume different competing configurations is the one with 
the lowest Gibbs free energy, G. For the present case, 
where the competition is among the DM and the AFM 
phases, the relevant quantity is the difference, AG = 
Gafm ~ Gdmj between their Gibbs energies. 



In general we find that the DM-AFM phase bound- 
ary moves in response to the disappearing of the dis- 
torted phase. Thus for the lower value of a/t (0.12), 
for which there is no distorted phase even at T = 0, the 
phase boundary does not change as the temperature is in- 
creased and Uc remains constant at ~ lAt. For the larger 
values of the el-ph interaction strengths investigated the 
T = phase diagram presents both a distorted and a 
homogeneous structural phase depending on U. In this 
case the DM-AFM phase boundary (i.e. Uc) decreases 
fast at low temperatures in response to the melting of the 
distorted phase and then becomes essentially constant. 

An important consequence of the movement of the 



AG = AH- TAS (6) 

where AH = Hafm — H^m and AS — Safm — Sum are 
respective the enthalpy and entropy differences. 

For standard spin crossover AH > 0, so that the most 
stable phase at low temperature in low spin (DM here). 
However, since the crossover transition is associated to 
the softening of the phonon modes of the first coordina- 
tion shell around the transition metal and to the forma- 
tion of a local magnetic moment, we also have AS' > 0. 
Hence as the temperature increases the entropic contri- 
bution to the Gibbs energy may eventually dominate over 
the enthalpic one and drive a phase transition. We now 
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want to establish that the same mechanism holds for the 
spin crossover region of the a-U diagram of Fig. [2j 

We have already demonstrated (see figure T]) that for 
T = the spin crossover region is occupied by the DM 
phase, meaning that AH > 0. Therefore one has only 
left to show that also AS* > 0. In general the entropy 
comprises of two main contributions, an electronic, 5"°', 
and a vibrational one, S^^^. Since the AFM phase is 
characterized by local spins, which are absent for the 
DM phase, we can immediately conclude that AS'"' — 
^AFM ~ '5'dm > 0. a more precise evaluation of S °' can 
be obtained by computing 



-fcBTr[p ln/5] 



(7) 



where p is the system density matrix. The calculated 
AS '^^ as a function of U jt for the representative value of 
a/t = 1, for which the spin crossover region is quite large, 
is presented in figure |4j The electron densities entering 
the evaluation of 5°' have been calculated as follows. For 
the low temperature phase (DM) p is calculated by fixing 
the occupation to Tij — 1/2 for every site and the geome- 
try is that obtained from the T = diagram. In contrast 
the entropy of the high-temperature phase is computed 
from a density matrix in which the occupation is fixed to 
the proper antiferromagnetic state (the temperature is 
T = O.OOli) and the geometry is again that of the T = 
solution. We have checked that the finite temperature 
geometry is rather similar to that obtained at T = for 
such a density-constraint solution. 

The general trend can be understood as follows. For 
small values of U jt there is no local magnetic moment 
formation, regardless of the chain geometry, so that 
AS"'' = 0. Then, as V jt gets larger a local magnetic 
moment gradually appears at each site, producing a lin- 
ear increase of AS'''. Such an increases continues until 
the local moment reaches the maximum value compati- 
ble with the chosen electron filling, which is 1 /ie (/ie is 
the Bohr magneton) for half-filling. At this point there 
is no further change in the electronic entropy and AS "' 
saturates to a positive value. 

Similarly we can also estimate the vibrational contri- 
bution to the entropy. This is obtained from the molecule 
phonon spectrum, fiwi, as 



S"*^ = fcR In 



1 



1 



+ fc] 



i/3 ^ huj^ 



are the same used for calculating the electronic contribu- 
tion to the entropy. Also in this case AS ^''^ as a function 
of U/t for a/t = 1 is presented in figure [4] 

In general the vibrational contribution to the entropy 
difference shows only a small dependence on the Coulomb 
on-site repulsion and approximately AS^''^ ^ 0. How- 
ever we also report a relatively sharp decreases as U/t 
approaches 1.7. This is value close to Uc for T = 0, i.e. 
corresponds to a region in the parameter space where our 
Monte Carlo analysis does not find a DM- AFM transi- 
tion and therefore system remains in the AFM state at 
any of the temperatures investigated. 

In summary the picture emerging from figure [4] is that 
of a region 0.84i < U < 1.85t in which AS is always 
positive. This is the only region of the parameter space 
where the entropy can drive the spin crossover and sub- 
stantially matches the spin crossover region observed in 
our Monte Carlo simulations (see Fig. [2]) . We then con- 
clude that also in this case, where the magnetic moment 
is not associated to the d shell of a transition metal, the 
spin crossover is entropy driven. 




(8) 



FIG. 4: (Color online) Entropy difference, AS, between the 
DM and AFM phases. The individual electronic, AS"', and 
vibrational, AS^''', contributions as well as their sum, AS, 
are presented as a function of U/t for a = It. The entropy 
is presented in adimensional units where both the energy and 
the temperature are in units of the hopping integral t. Note 
that there is a large region 0.84t < U < 1.85t, where AS > 
0. In this region the entropic contribution to the Gibbs free 
energy can drive a spin crossover transition. 



where Wj is the vibration frequency of the i-th mode 
and Njj is the total number of modes. We calculate the 
phonon spectra of the AFM and DM configurations by 
diagonalizing the associated dynamical matrices. These 
are constructed by finite difference, i.e. by displacing the 
atomic sites by a small fraction of the equilibrium bond 
length (0.1 %) and then by numerically evaluating the 
energy gradient (the force) associated to such a displace- 
ment. The density matrices and the initial geometries 
used to construct the finite difference dynamical matrices 



Finally before summarizing, we wish to comment on 
the role played by the hyperfine interaction. In general 
we expect very little changes to the phase diagram ob- 
tained by neglecting the hyperfine contribution to the 
total energy, since this is rather small for realistic hy- 
perfine coupling strengths. In particular we have verified 
that minor modifications to the T = phase diagram 
start to appear for Jh is the region of 0.01 1, which cor- 
responds to local magnetic fields of 10^ T (considering 
I Si I = 1 and t of the order of 1 eV) . 
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IV. CONCLUSION 

In summary we have discussed the phase diagram of 
an extended Hubbard-Peierls model chosen to emulate 
the electronic structure of tt— conjugated polymer chains. 
The model has been explored both in its T = ground 
state and at finite temperature by Monte Carlo methods 
combined with a mean field treatment of the Hubbard 
many-body interaction. At T = the model presents 
four different phases depending on the relative strength 
of the Coulomb on-site repulsion U/t and the el-ph cou- 
pling strength a/t. The four phases are characterized ei- 
ther by a diamagnetic or an antiferromagnetic magnetic 
state and by the possible presence of a dimerized geo- 
metrical configuration. By increasing the temperature 
the structural distortion disappears and only the phase 
boundary between the DM and the AFM solution re- 
mains. Intriguingly the position in the a-U parameter 
space of such a phase boundary changes with temper- 



ature so that there is a region of parameters where a 
temperature driven DM-AFM spin crossover transition 
can be found. We have investigated the nature of the 
phase transition by calculating the relative entropy of 
the different magnetic phases and found that this is in- 
deed entropy driven as in the most conventional case of 
divalent magnetic molecules. This suggests that 7r-type 
magnetism can be achieved in organic polymers and that 
this can be tuned by temperature. 
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